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1 Introduction 

In this paper we present a fast algorithm for the numerical solution of systems 
of reaction-diffusion equations, 

dtu + a ■ Vit = Au + F(x,t,u), x G n c R 3 , t > 0. (1) 

Here, u is a vector-valued function, u = u(x, i) G R m , m is large, and the 
corresponding system of ODEs, dtu = F(x,t,u), is stiff. Typical examples 
arise in air pollution studies, where a is the given wind field and the nonlinear 
function F models the atmospheric chemistry. 

The time integration of Eq. (|l|) is best handled by the method of character- 
istics pL The problem is thus reduced to designing for the reaction-diffusion 
part a fast solver that has good stability properties for the given time step and 
does not require the computation of the full Jacobi matrix. 

An operator-splitting technique, even a high-order one, combining a fast 
nonlinear ODE solver with an efficient solver for the diffusion operator is less 
effective when the reaction term is stiff. In fact, the classical Strang splitting 
method may underperform a first-order source splitting method [Q. The al- 
gorithm we propose in this paper uses an a posteriori filtering technique to 
stabilize the computation of the diffusion term. The algorithm parallelizes well, 
because the solution of the large system of ODEs is done pointwise; however, the 
integration of the chemistry may lead to load-balancing problems [|[ ||. The 
Tchebycheff acceleration technique proposed in 0] offers an alternative that 
complements the approach presented here. 

To facilitate the presentation, we limit the discusssion to domains that ei- 
ther admit a regular discretization grid or decompose into subdomains that ad- 
mit regular discretization grids. We describe the algorithm for one-dimensional 
domains in Section 2 and for multidimensional domains in Section 3. Section 4 
briefly outlines future work. 
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2 One-Dimensional Domains 



Consider the scalar equation 

d t u = d 2 x u + f(u), x G (0,7r), t > 0. 



(2) 



We combine a backward Euler approximation in time with an explicit finite- 
difference approximation of the diffusive term, 

— = 2D xx u n - D x u n - X + f{u n+1 ). (3) 
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This scheme is second-order accurate in both space and time |], ^[ [| . To analyze 
its stability, we take the Fourier transform of the linear equation, 
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where = 2h 2 (cos(/ifc) — 1), from which we obtain the stability condition 



2^ 

K 2 





- 1 


cos 


\N ) 





4 , 7T 

< -, h = — . 

3' N 



Thus we conclude that the time step must satisfy the constraint 

At < \h 2 . 



(5) 
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However, this constraint is imposed by the high frequencies, which are poorly 
handled by second-order finite differences anyway. For example, with central 
differences, the relative error for high-frequency waves cos(fcx) with k w N can 
grow at a rate of up to 9%. The idea is therefore to relax the constraint on 
the time step by applying a filter after each time step, which removes the high 
frequencies but maintains second-order accuracy in space. 



2.1 Filters 

By a filter of order p we mean an even function a : R — > R that satisfies the 
conditions (i) ct(0) = 1, (ii) <jW(0) = for I = 1, . . . ,p - 1, (iii) a(rj) = for 
M > 1, and (iv) a <E CP-^R). 

Theorem M . Let / be a piecewise C p function with one point of discontinuity, 
£, and let a be a filter of order p. For any point y G [0, 2tt], let d(y) = min{|y — 
£ + 2kn\:k= -1,0,1}. If f a N -Efc-oo fk<r(k/N)e ik v, thcn 

|/(y) - k/N\ < CN^idiy^-VKif) + CN^-Py^Us, 

where 

P- 1 poo 

K(f) = - / (0 (Di / iGf-'^jidf,. 

In other words, a discontinuity of / leads to a Fourier expansion with an error 
that is O(l) near the discontinuity and 0(N~ 1 ) away from the discontinuity. 
We must therefore apply a shift and extend to [0, 2ir] before applying a filter. 
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2.2 The Algorithm 

We now describe the postprocessing algorithm that is to be applied after each 
time step. (We do not explicitly indicate the dependence of u on the time step, 
and we use the abbreviations u Q — it(0) and — u(tt).) 
First, we apply a low-frequency shift, 

v(x) = u(x) — (ai + «2 cos(x)), ai — ^(u — u„), a 2 — ^(u + u n ). (7) 

Then we extend v to (0, 27r), using the definition 

v(2ir -x) = -v(x), x € (0, tt). (8) 

Thus, v is a 27r-periodic function in C 1 (0,27r). Let v k be the fcth coefficient of 
its Fourier expansion. 

Next, we apply an eighth-order filter ||, 

„ L M A „ikx 



a N v{x)=^aU—\v k e lkx , (9) 

k ^ ' 

where 

<r(0 - (35 - 84y + 70y 2 - 20y 3 )y 4 , y = 2/(0 = |(1 + cos(ttO). (10) 

Here, k is a stretching factor, k > 1. The correct choice of k follows from a 
Fourier analysis of Eq. (||) , 

7T 

K>Kc = cos(l-2/ 1 2 /(3A0)' 

The choice k — ^k c gives satisfactory results, but in principle one can compute 
the optimum value of k at each time step by monitoring the growth of the 
high-frequency waves that have not been completely filtered out. 
Finally, we recover u from the inverse shift, 

u(x) = <tnv{x) + a\ + Q!2 cos(a;). (12) 

The theorem quoted in the preceding section shows that the filtering process 
may affect the spatial accuracy of the method. Since the filter is applied to a 
27r-periodic function that is C 1 at the points Xk = kir, k € Z, and C 2 everywhere 
else, the error is of the order of N~ 2 in the neighborhood of x k and iV~ 3 away 
from Xk- In principle, we maintain therefore second-order accuracy in space as 
long as k is of order one. 

If the solution is in C 3 (0, it) at each time level, we can improve the algorithm 
by replacing the first-order shift (|7|) by a third-order shift, 

4 

v(x) — u(x) — ^ aj cos((j — l)x), (13) 



3 



such that the extension of v to a 27r-periodic function is in C 3 (0, 2tt). The first- 
and third-order derivatives of v are zero at the points x k , and the second-order 
derivative is approximately given by 
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f(u n+1 (x k )). 



(14) 



The coefficients aj are found by solving a linear system of equations, 

Q'l + «2 + Oi 3 + Q!4 = u(0), 

a'l - a 2 + a 3 - a 4 — u(ir), 
-a 2 - 4a 3 - 9a 4 = u xx (0), 
a 2 - 4a 3 + 9a 4 = u xs (7r). 

The third-order shift improves the performance of the filter for large k and 
allows for a larger time step. 

2.3 Numerical Results 

Figure [l] shows some accuracy results for Eq. (||) , where 

u(x,t) = cos(t)({x/ir) 4 +cos(3a;)), x e (0, tt), t > 0. (15) 
We observe a plateau for small time steps, when the second-order spatial er- 
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Figure 1: Accuracy of the stabilized explicit scheme for the heat equation. The 
horizontal coordinate is 3Ai//i 2 . * : first-order shift; o : third-order shift. 



ror dominates. The second-order error in time becomes dominant as the time 
step increases. The figure confirms the superior performance of the third-order 
shift (|l^) over the first-order shift (^) at large time steps. 
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Although the algorithm is based only on linear stability considerations, it is 
still effective for systems of nonlinear reaction-diffusion equations. In Figure |^ 
we present some results for a predator-prey system, 



dtu = d xx u + au- buv, dtv = d xx v — cu — duv, x G (0, tt), t > 0, (16) 
with a = 1.2, b = 1.0, c = 0.1, and d = 0.2. At these parameter values, the 
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Figure 2: Accuracy of the stabilized explicit scheme for the predator-prey sys- 
tem. The horizontal coordinate is 3At/h 2 . * : first-order shift; o : third-order 
shift. 



ODE system (reactions only) has a limit cycle. However, when the boundary 
conditions are constant in time, the solution of the reaction-diffusion system 
goes to steady state. To build a relevant test case for the algorithm, we impose 
periodic excitations at both boundaries, 

u(0/-7r, t) = Uo/7r(l + cos(i)), v(0/w,t) = t>o/7r(l + cos(i)), t > 0. 

Although the time step can still be limited by nonlinear instabilities, we never 
observed negative values of the unknowns u and v, which are commonly associ- 
ated with such instabilities. 

The algorithm (Q)— ( 12 ) extends in a straightforward way when one uses a 



domain-decomposition scheme with overlapping subdomains. One simply ap- 
plies the same algorithm at each time step to each subdomain separately. How- 
ever, the number of waves per subdomain is of the order of the total number 
of grid points, N, divided by the number of subdomains, Nd, so the balance 
between the order of accuracy of the filter — (N/(kN,i))~ 3 for a first-order filter 
or (N/ (kN^))^ 5 for a third-order filter — and the second-order accuracy N~ 2 
of the spatial discretization of the underlying algorithm (||) deteriorates as Nd 
increases. The maximum time step for which the scheme remains stable may 
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become even less than when no domain decomposition is used. Furthermore, the 
Gibbs phenomenon tends to destabilize the algorithm. This phenomenon is a 
consequence of the jump in the derivatives at the endpoints of the subdomains 
(second-order derivatives in the case of the first-order shift ([?]), fourth-order 
derivatives in the case of the third-order shift (|l3|)). Since the Gibbs phe- 
nomenon arises at the artificial interface and is damped away from it, an increase 
of the overlap generally produces a composite signal u that has fewer oscillations 
than each of the piecewise (overlapping) components. One can therefore obtain 
good results by adapting the size of the overlap. The larger the overlap, the 
larger the time step that can be taken; see Figure 0. 



* for d=1 and N=63, o for d=3 and N=61 , + for d=5 and N=63 
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Figure 3: Accuracy of the stabilized explicit scheme applied to the heat equation 
with overlapping subdomains. 



3 Two-Dimensional Domains 

We now consider a Dirichlet problem in two dimensions, 

d t u = Au + f(u), (x,y) e (0,tt) 2 , t >0, (17) 
u{x,0/tt) = g 0/n (x), u(0/n,y) = h 0/7I (y), x, y € (0, tt), (18) 

where the functions g satisfy the compatibility conditions <?o/7r(0) = ho(0/n) 
and ffo/7r( 7r ) = h^iQ/ir). We consider a numerical scheme similar to Eq. (||), 
where the diffusive term is approximated, for example, by a five-point stencil. 
The postprocessing algorithm is essentially the same, except that we need an 
apropriate low-frequency shift so we can apply a filter to a smooth periodic 
function in two space dimensions. The shift is constructed in two steps. In the 
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first step, we render the boundary condition in the x direction homogeneous, 



v(x,y) = u(x,y) - {a x {y) + a 2 (y) cos(x)), 

ai{y) = Kfo(y) -9n{y)), a 2 (y) = \[go{y) + 9^{v))- 



(19) 
(20) 



In the second step, we shift in the y direction, 



w(x,y) = v(x,y) - ((3i(x) + (3 2 {x) cos(y)), 
/3i(x) = |(«(a;,0) -v(x,n)), (3 2 {x) = ±(v(x,0) + v(x,w)). 



(21) 
(22) 



The final step is the reconstruction step, 

u(x,y) = a N w(x,y) + ai(y) + a 2 (y)cos(x) + j3i{x) + 2 {x) cos(y). (23) 

To make sure that no high-frequency waves remain, we filter the high-frequency 
components from the boundary conditions g with a procedure similar to ((7])- 



It is much more difficult to construct a high-order filter similar to (|13J) in 
two dimensions, because the second-order derivatives cannot be obtained from 
the PDE, as in the one-dimensional case (|l4|). So far, we have used only the 
first-order shifts (|l9|) and (|2l]) in our numerical experiments. Nevertheless, the 
algorithm allows for a significant increase of the time step. We have also tested 
the domain-decomposition version of the algorithm, using strip subdomains with 
an adaptive overlap, with good results. 

We note that the computation in each block can be done in parallel and 
that the Jacobi matrix does not depend on the spatial variables. The arithmetic 
complexity of the algorithm is therefore relatively small. Also, the algorithm is 
suitable for multicluster architectures. Each block can be assigned to a cluster, 
and parallel fast sine transforms can be used for the filtering process inside each 
cluster. The cost of communication between blocks is minimal, since the scheme 
is similar to the communication scheme of the additive Schwarz algorithm. 

4 Conclusion 

In this paper we have presented a postprocessing algorithm that stabilizes the 
time integration of systems of reaction-diffusion equations when the diffusion 
term is treated explicitly. The algorithm is easy to code and can be combined 
with domain-decomposition methods that use regular grids in each subblock. 
In future work, we will consider the performance of its parallel implementation 
and its robustness for large systems of reaction-diffusion equations with stiff 
chemistry, which arise in some air pollution models. 
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